import pandas as pd
import numpy as np
import seaborn as sns
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.graphics.regressionplots import plot_partregress_grid
import matplotlib.pyplot as plt
%matplotlib inline

Core Analysis¶

After data cleaning, we visualized our data with the matrix plots first.

systolic blood pressure¶

pair_plot_data_sys_bp = full_dt[['avg_BPXSY', 'gender', 'age', 'bmi', 'sleep', 
                                 'smoke', 'workhrs', 'avg_alcohol_freq_wk']]
ax = pd.plotting.scatter_matrix(pair_plot_data_sys_bp, alpha = 0.5, figsize = (14, 12)); 
plt.xticks(fontsize = 2)
plt.yticks(fontsize = 2);

diatolic blood pressure¶

pair_plot_data_dia_bp = full_dt[['avg_BPXDI', 'gender', 'age', 'bmi', 'sleep', 
                                 'smoke', 'workhrs', 'avg_alcohol_freq_wk']]
pd.plotting.scatter_matrix(pair_plot_data_dia_bp, alpha = 0.5, figsize = (14, 12)) 
plt.xticks(fontsize = 2)
plt.yticks(fontsize = 2);

As shown in the matrix plots for the systolic blood pressure, we found that there were weak correlations between the variables. As for diastolic blood pressure, the matrix plots show that there are weak relations or almost no relations between the variables. Moreover, from the histograms, we found that some variables were not normally distributed; however, most variables, including systolic and diastolic blood pressure, were normally distributed.

After the preliminary analysis on the data, we fitted our linear regression models.

systolic blood pressure¶

For the systolic blood pressure, we fit the model: $$ \begin{aligned} y &= \beta_{0} + \beta_{workhrs}*x_{workhrs} + \beta_{gender}*x_{gender} + \beta_{age}*x_{age} + \beta_{bmi}*x_{bmi} \notag\\ &+ \beta_{sleep}*x_{sleep} + \beta_{smoke}*x_{smoke} + \beta_{avg\_alcohol\_freq\_wk}*x_{avg\_alcohol\_freq\_wk} + \epsilon \end{aligned} $$

model_formula_sys_bp = 'avg_BPXSY ~ gender + age + bmi + sleep + smoke + workhrs + avg_alcohol_freq_wk'
model_sys_bp = smf.ols(formula = model_formula_sys_bp, data = full_dt)
model_fit_sys_bp = model_sys_bp.fit()

print(model_fit_sys_bp.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:              avg_BPXSY   R-squared:                       0.089
Model:                            OLS   Adj. R-squared:                  0.061
Method:                 Least Squares   F-statistic:                     3.192
Date:                Mon, 09 Dec 2019   Prob (F-statistic):            0.00303
Time:                        18:28:00   Log-Likelihood:                -1037.0
No. Observations:                 238   AIC:                             2090.
Df Residuals:                     230   BIC:                             2118.
Df Model:                           7                                         
Covariance Type:            nonrobust                                         
=======================================================================================
                          coef    std err          t      P>|t|      [0.025      0.975]
---------------------------------------------------------------------------------------
Intercept              90.2929     11.828      7.634      0.000      66.988     113.598
gender                 -3.6250      2.744     -1.321      0.188      -9.031       1.781
age                     0.4704      0.127      3.697      0.000       0.220       0.721
bmi                     0.6472      0.202      3.196      0.002       0.248       1.046
sleep                   0.4919      0.980      0.502      0.616      -1.439       2.423
smoke                   1.1942      2.651      0.450      0.653      -4.029       6.418
workhrs                 1.5753      2.677      0.588      0.557      -3.699       6.850
avg_alcohol_freq_wk     0.0192      0.720      0.027      0.979      -1.399       1.437
==============================================================================
Omnibus:                       35.325   Durbin-Watson:                   1.945
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               52.892
Skew:                           0.877   Prob(JB):                     3.27e-12
Kurtosis:                       4.502   Cond. No.                         522.
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

diatolic blood pressure¶

For the systolic blood pressure, we fitted the following model: $$ \begin{aligned} y &= \beta_{0} + \beta_{workhrs}*x_{workhrs} + \beta_{gender}*x_{gender} + \beta_{age}*x_{age} + \beta_{bmi}*x_{bmi} \notag\\ &+ \beta_{sleep}*x_{sleep} + \beta_{smoke}*x_{smoke} + \beta_{avg\_alcohol\_freq\_wk}*x_{avg\_alcohol\_freq\_wk} + \epsilon \end{aligned} $$

model_formula_dia_bp = 'avg_BPXDI ~ gender + age + bmi + sleep + smoke + workhrs + avg_alcohol_freq_wk'
model_dia_bp = smf.ols(formula = model_formula_dia_bp, data = full_dt)
model_fit_dia_bp = model_dia_bp.fit()

print(model_fit_dia_bp.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:              avg_BPXDI   R-squared:                       0.092
Model:                            OLS   Adj. R-squared:                  0.065
Method:                 Least Squares   F-statistic:                     3.347
Date:                Mon, 09 Dec 2019   Prob (F-statistic):            0.00204
Time:                        18:28:00   Log-Likelihood:                -954.08
No. Observations:                 238   AIC:                             1924.
Df Residuals:                     230   BIC:                             1952.
Df Model:                           7                                         
Covariance Type:            nonrobust                                         
=======================================================================================
                          coef    std err          t      P>|t|      [0.025      0.975]
---------------------------------------------------------------------------------------
Intercept              49.4133      8.348      5.919      0.000      32.965      65.862
gender                 -4.5041      1.937     -2.326      0.021      -8.320      -0.688
age                     0.2760      0.090      3.073      0.002       0.099       0.453
bmi                     0.5378      0.143      3.763      0.000       0.256       0.819
sleep                   0.3394      0.692      0.491      0.624      -1.024       1.702
smoke                   0.3443      1.871      0.184      0.854      -3.342       4.031
workhrs                 0.2246      1.889      0.119      0.905      -3.498       3.947
avg_alcohol_freq_wk     0.0909      0.508      0.179      0.858      -0.910       1.092
==============================================================================
Omnibus:                       12.349   Durbin-Watson:                   2.163
Prob(Omnibus):                  0.002   Jarque-Bera (JB):               17.669
Skew:                          -0.352   Prob(JB):                     0.000146
Kurtosis:                       4.134   Cond. No.                         522.
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

From the summary results of the linear regression model for the systolic blood pressure, we found the coefficient of the working hours is positive, which implies that working overtime increases systolic blood pressure level; however, this effect is not significant ($\beta = 1.58$, p = 0.557). As a result, we can conclude that working overtime is not significantly correlated to systolic blood pressure. Similarly, from the summary results of the linear regression model for diastolic blood pressure, we found that the coefficient of the working hours is positive, which implies that working overtime increases diastolic blood pressure level; however, this effect was not significant as well ($\beta = 0.2246$, p = 0.905). As a result, we conclude that working overtime is not significantly correlated to diastolic blood pressure.

After fitting the linear regression model, we checked whether the following assumptions of the linear regression model are held:

Homoscedasticity: The variance of the residual is constant.
Linearity: The relationship between the independent and dependent variables is linear.
No or little multicollinearity: There is no or little collinearity between independent variables.

In order to check the homoscedasticity assumption, we plotted standardized residuals against fitted values plots for both systolic blood pressure and diastolic blood pressure.

systolic blood pressure¶

# fitted values (need a constant term for intercept)
model_fitted_y_sys_bp = model_fit_sys_bp.fittedvalues

# model residuals
model_resid_sys_bp = model_fit_sys_bp.resid

# standardized residuals
model_std_resid_sys_bp = model_fit_sys_bp.get_influence().resid_studentized_internal

fig, ax = plt.subplots()
sns.residplot(model_fitted_y_sys_bp, model_std_resid_sys_bp, lowess = True,
              scatter_kws = {'alpha': 0.5}, 
              line_kws = {'color': 'red', 'lw': 1, 'alpha': 0.8})
plt.title('standardized residual plot (systolic blood pressure)')
ax.set_xlabel('fitted values')
ax.set_ylabel('standardized residuals')

<matplotlib.text.Text at 0x112b84550>

diatolic blood pressure¶

# fitted values (need a constant term for intercept)
model_fitted_y_dia_bp = model_fit_dia_bp.fittedvalues

# model residuals
model_resid_dia_bp = model_fit_dia_bp.resid

# standardized residuals
model_std_resid_dia_bp = model_fit_dia_bp.get_influence().resid_studentized_internal

# sqrt absolute standardized residuals
model_abs_sqrt_std_resid_dia_bp = np.sqrt(np.abs(model_std_resid_dia_bp))

# absolute residuals
model_abs_resid_dia_bp = np.abs(model_resid_dia_bp)

fig, ax = plt.subplots()
sns.residplot(model_fitted_y_dia_bp, model_std_resid_dia_bp, lowess = True,
              scatter_kws = {'alpha': 0.5}, 
              line_kws = {'color': 'red', 'lw': 1, 'alpha': 0.8})
plt.title('standardized residual plot (diatolic blood pressure)')
ax.set_xlabel('fitted values')
ax.set_ylabel('standardized residuals')

<matplotlib.text.Text at 0x1128f4668>

From the standardized residuals plot of the systolic blood pressure, we could find the mean of the residual was roughly 0, and the variance was roughly constant. Also, we had the same findings for diastolic blood pressure. Therefore, the assumption of homoscedasticity could be considered as satisfied.

In order to check the linearity assumption, we plotted partial regression plots for both systolic blood pressure and diastolic pressure against each of the independent variables.

systolic blood pressure¶

fig = plt.figure(figsize=(8, 12))
plot_partregress_grid(model_fit_sys_bp, fig=fig)
plt.show()

diatolic blood pressure¶

fig = plt.figure(figsize=(8, 12))
plot_partregress_grid(model_fit_dia_bp, fig=fig)
plt.show()

From the partial regression plots of the systolic blood pressure, we could find that for each of the independent variables, the expected value of the dependent variable (systolic blood pressure) was indeed a linear function of the independent variable, controlled by all the other variables. We can get the same conclusion from the partial regression plots of diastolic blood pressure. Therefore, the assumption of linearity can be considered as satisfied.

In order to check the multicollinearity, we computed the Pearson correlations between each two continuous independent variables and the Pearson correlations between each two binary independent variables.

check correlations between the continuous variables¶

continuous_vars = full_dt[['age', 'bmi', 'sleep', 'avg_alcohol_freq_wk']]
corr_continuous = np.corrcoef(continuous_vars.values.T)
# corr_continuous

fig, ax = plt.subplots(figsize=(12, 10))
sns.heatmap(corr_continuous, annot=True, square=True, annot_kws = {"size": 12})
plt.title("correlations between continuous variables", fontsize = 16)
ax.set_xticklabels(continuous_vars)
ax.set_yticklabels(continuous_vars, va = 'center')
ax.tick_params(axis='x', labelsize=12)
ax.tick_params(axis='y', labelsize=12)

check correlations between categorical variables¶

categorical_vars = full_dt[['gender', 'smoke', 'workhrs']]
categorical_vars.apply(lambda x : pd.factorize(x)[0]).corr(method='pearson')

As shown in the heatmap of the correlations between the continuous independent variables, we could find that there are very small or almost no correlations between the continuous independent variables. Also, as the table of the Pearson chi-squared correlation test results showed, we could find that there is little collinearity between each pair of the binary variables.

Additional Analysis¶

Here, I fit the logistic regression (as additional analysis) to examine the association between the predictors and the risk of having the abnormal blood pressure level. Following the guide from the American Heart Association, the abnormal systolic and diastolic blood pressure were coded as 1 if the respondent had systolic blood pressure higher than 140mm Hg or diastolic blood pressure higher than 90mm Hg. If not, they were coded 0 to indicate normal blood pressure level.

# convert the variables systolic blood pressure and diatolic blood pressure into binary variables
full_dt['bin_avg_BPXSY'] = np.where(full_dt.avg_BPXSY >= 140, 1, 0)
full_dt['bin_avg_BPXDI'] = np.where(full_dt.avg_BPXDI >= 90, 1, 0)

# drop the variables avg_BPXSY and avg_BPXDI which are not needed for the logistic models
add_analysis_dt = full_dt.drop(['id', 'race', 'avg_BPXSY', 'avg_BPXDI'], axis = 1)
add_analysis_dt.head()

systolic blood pressure¶

endog_sys_bp = add_analysis_dt.bin_avg_BPXSY
exog_sys_bp = sm.add_constant(add_analysis_dt.iloc[:, 0:7])
logit_model_sys_bp = sm.Logit(endog_sys_bp, exog_sys_bp)
logit_model_fit_sys_bp = logit_model_sys_bp.fit()
print(logit_model_fit_sys_bp.summary())

Optimization terminated successfully.
         Current function value: 0.589699
         Iterations 6
                           Logit Regression Results                           
==============================================================================
Dep. Variable:          bin_avg_BPXSY   No. Observations:                  238
Model:                          Logit   Df Residuals:                      230
Method:                           MLE   Df Model:                            7
Date:                Mon, 09 Dec 2019   Pseudo R-squ.:                 0.06776
Time:                        18:28:04   Log-Likelihood:                -140.35
converged:                       True   LL-Null:                       -150.55
                                        LLR p-value:                  0.004764
=======================================================================================
                          coef    std err          z      P>|z|      [0.025      0.975]
---------------------------------------------------------------------------------------
const                  -6.2485      1.500     -4.165      0.000      -9.189      -3.308
gender                 -0.1850      0.318     -0.581      0.561      -0.809       0.439
age                     0.0531      0.016      3.343      0.001       0.022       0.084
bmi                     0.0684      0.024      2.872      0.004       0.022       0.115
sleep                   0.1249      0.116      1.078      0.281      -0.102       0.352
smoke                   0.0841      0.307      0.274      0.784      -0.518       0.686
workhrs                 0.2389      0.304      0.785      0.432      -0.357       0.835
avg_alcohol_freq_wk     0.0682      0.081      0.846      0.397      -0.090       0.226
=======================================================================================

For the systolic blood pressure, the logistic regression model indicated that working overtime is related to increased risk of having abnomral systolic blood pressure ($\beta = 0.24$, p = 0.432); however, this association was not statistically significant. Among the covariates, higher age and higher bmi were significantly associated with increased risk of having abnormal systolic blood pressure.

diatolic blood pressure¶

endog_dia_bp = add_analysis_dt.bin_avg_BPXDI
exog_dia_bp = sm.add_constant(add_analysis_dt.iloc[:, 0:7])
logit_model_dia_bp = sm.Logit(endog_dia_bp, exog_dia_bp)
logit_model_fit_dia_bp = logit_model_dia_bp.fit()
print(logit_model_fit_dia_bp.summary())

Optimization terminated successfully.
         Current function value: 0.414587
         Iterations 6
                           Logit Regression Results                           
==============================================================================
Dep. Variable:          bin_avg_BPXDI   No. Observations:                  238
Model:                          Logit   Df Residuals:                      230
Method:                           MLE   Df Model:                            7
Date:                Mon, 09 Dec 2019   Pseudo R-squ.:                 0.08441
Time:                        18:28:04   Log-Likelihood:                -98.672
converged:                       True   LL-Null:                       -107.77
                                        LLR p-value:                   0.01113
=======================================================================================
                          coef    std err          z      P>|z|      [0.025      0.975]
---------------------------------------------------------------------------------------
const                  -7.0335      1.873     -3.754      0.000     -10.705      -3.362
gender                 -0.6208      0.419     -1.482      0.138      -1.442       0.200
age                     0.0522      0.020      2.576      0.010       0.012       0.092
bmi                     0.0862      0.029      2.931      0.003       0.029       0.144
sleep                   0.0128      0.144      0.089      0.929      -0.269       0.295
smoke                   0.2854      0.383      0.746      0.456      -0.465       1.036
workhrs                 0.4136      0.375      1.104      0.270      -0.321       1.148
avg_alcohol_freq_wk     0.1144      0.094      1.211      0.226      -0.071       0.299
=======================================================================================

For the diastolic blood pressure, the logistic regression model showed that working overtime is related to increased risk of having abnomral diastolic blood pressure ($\beta = 0.41$, p = 0.27); however, this association was not statistically significant. Among the covariates, higher age and higher bmi were significantly associated with increased risk of having abnormal diastolic blood pressure.

	gender	smoke	workhrs
gender	1.000000	0.123827	-0.077188
smoke	0.123827	1.000000	0.120028
workhrs	-0.077188	0.120028	1.000000

	gender	age	bmi	sleep	smoke	workhrs	avg_alcohol_freq_wk	bin_avg_BPXSY	bin_avg_BPXDI
0	0	46.0	27.6	10.0	1	0	6.000000	1	1
1	1	54.0	38.8	8.0	0	0	0.920598	0	0
2	1	45.0	33.2	8.0	0	1	6.000000	0	1
3	0	53.0	23.3	7.0	1	0	2.000000	1	0
4	0	52.0	27.2	8.0	0	0	1.000000	1	1

Project_Group5

Hyesue Jang, Tian Wang, Zhe Zhao

Dec 11, 2019

Introduction

Data Selction and Cleaning

Method

Core Analysis and Results

lm and R_data.table

additional analysis

GLS and R_dplyr

Additional Analysis (Hyesue Jang)

Python

Core Analysis¶

systolic blood pressure¶

diatolic blood pressure¶

systolic blood pressure¶

diatolic blood pressure¶

systolic blood pressure¶

diatolic blood pressure¶

systolic blood pressure¶

diatolic blood pressure¶

check correlations between the continuous variables¶

check correlations between categorical variables¶

Additional Analysis¶

systolic blood pressure¶

diatolic blood pressure¶

References